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array  response  an  improved  estimate  of  the  incident  wave- 
number  spectrum  may  be  obtained.  In  particular  the  method  of 
progressive  substitutions  has  been  widely  used  to  deconvolve 
the  array  window.  By  assuming  all  the  incident  signals  to  be 
a  linear  superposition  of  a  finite  number  of  plane  waves  this 
method  is  extended  to  deconvolve  the  narrowband  beam  outputs  of 
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iterative  solutions  is  proved  and  the  relationship  of  the  limits 
to  some  recently  proposed  quadratic  estimators  is  shown. 


Approved  for  Public  Release 


POSTAL  ADDRESS:  Chief  Superintendent,  Weapons  Systems  Research  Laboratory, 

Box  2151,  G.P.O.,  Adelaide,  South  Australia,  5001. 


UNCLASSIFIED 


DOCUMENT  CONTROL  DATA  SHEET 


Security  classification  of  this  page 


DOCUMENT  NUMBERS 


UNCLASSIFIED 

2  SECURITY  CLASSIFICATION 


AR 

Number: 

AR-001-954 

Report 

Number: 

WSRL-0141-TR 

V 

Other 

Numbers: 

a.  Complete 
Document : 

Unclassified 

b.  Title  in 

Isolation: 

Unclassified 

c.  Summary  in 

Isolation: 

Unclassified 

3  TITLE 


AN  ITERATIVE  APPROACH  TO  THE  DECONVOLUTION  OF  THE 
NARROWBAND  OUTPUTS  OF  A  CONVENTIONAL  BEAMFORMER 


5  DOCUMENT  DATE: 


February  1980 


6 

6.1 

TOTAL  NUMBER 
OF  PACES 

22 

6.2 

NUMBER OF 
REFERENCES: 

6 

7  7.1  CORPORATE  AUTHOR(S): 


Weapons  Systems  Research  Laboratory  v- 


7.2  DOCUMENT  SERIES 
AND  NUMBER 

Weapons  Systems  Research  Laboratory 
0141-TR 


8  REFERENCE  NUMBERS 


II  I  COMPUTER  PROGRAM(S) 
(Title(s)  and  language(s)) 


10  IMPRINT  (Publishing  organisation) 


Defence  Research  Centre  Salisbury 


1 2  RELEASE  LIMITATIONS  (of  the  document): 


Approved  for  Public  Release 


1 2.0  OVERSEAS 


UNCLASSIFIED 


Security  classification  of  this  page:  I  UNCLASSIFIED 


14 

a. 


DESCRIPTORS: 

EJC  Thesaurus 
Terms 


Antenna  radiation  patterns 
Adaptive  systems 
Antenna  arrays 


b.  Non-Thesaurus 
Terms 


Beamforming  Adaptive  techniques 

Polar  response  Pseudo  inverse 

Deconvolution 
Estimation  theory 


15 


COSATI  CODES 


2014 


16 


LIBRARY  LOCATION  CODES  (for  libraries  listed  in  the  distribution): 


17  |  SUMMARY  OR  ABSTRACT: 

(if  this  is  security  classified,  the  announcement  of  this  report  will  be  similarly  classified) 


In  conventional  beamforming  the  estimated  noise  power 
wave  number  spectrum  is  the  convolution  of  the  incident 
distribution  with  the  array  response.  By  removing  the 
effect  of  the  array  response  an  improved  estimate  of  the 
incident  wave-number  spectrum  may  be  obtained.  In 
particular  the  method  of  progressive  substitutions  has 
been  widely  used  to  deconvolve  the  array  window.  By 
assuming  all  the  incident  signals  to  be  a  linear  super¬ 
position  of  a  finite  number  ‘of  plane  waves  this  method  is 
extended  to  deconvolve  the  narrowband  beam  outputs  of  a 
conventional  beamformer.  Two  methods  are  used,  the  first 
to  deconvolve  the  complex  beam  amplitudes  and  the  second 
to  deconvolve  the  beam  powers.  The  convergence  of  these 
proposed  iterative  solutions  is  proved  and  the  relationship 
of  the  limits  to  some  recently  proposed  quadratic  estimators 
is  shown. 


Security  classification  of  this  page: 


UNCLASSIFIED 


WSRL-0141 


TABLE  OF  CONTENTS 

Page  No 


1.  INTRODUCTION  1  -  2 

2.  FORMULATION  OF  THE  ITERATIVE  METHOD  2-6 

2.1  Deconvolution  of  complex  beam  outputs  3-5 

2.2  Deconvolution  of  beam  powers  5-6 

3.  SOLUTION  AND  CONVERGENCE  OF  THE  ITERATIVE  METHOD  7-9 

3.1  Convergence  7 

3.2  Limit  (complex  beam  outputs)  7-8 

3.3  Limit  (beam  powers)  9 

4.  DISCUSSION  10-13 

4.1  Complex  beam  outputs  10  -  12 

4.2  Beam  powers  12  -  13 

5.  EXAMPLES  13-16 

5.1  K  receivers  and  1  plane  wave  15 

5.2  K  receivers  and  N  sources  at  half  wavelength  spacing  15 

5.3  K  receivers  and  K  sources  at  quarter  wavelength  spacing  15  -  16 

6.  CONCLUSIONS  16 

REFERENCES  17 

LIST  OF  APPENDICES 

I  A  VECTOR  IDENTITY  18-19 

II  CONVERGENCE  OF  A  SERIES  20  -  21 


Figure  1.  N  beams  steered  in  the  physical  region  at  d  =  j 


A 


Di  j: 


WSRL-0141-TR 


-  1  - 


1 .  INTRODUCTION 

In  conventional  beamforming  it  is  well  known  that  the  estimated  complex  wave- 
number  spectrum  of  the  incident  noise  distribution  is  the  convolution  of  the 
true  wave-number  spectrum  with  the  Fourier  transform  of  the  array  window. 

Since  the  array  window  is  zero  outside  a  finite  region,  a  consequence  of  this 
convolution  is  to  limit  the  wave-number  (and  hence  angular)  resolution. 
Deconvolution  techniques  can  be  used  to  remove  (or  minimise)  the  effect  of  the 
array  window  and  consequently  improve  the  estimate  of  the  wave-number  spectrum. 
Unfortunately  deconvolution  methods  do  not  always  provide  a  unique  estimate  of 
the  wave-number  spectrum  and  constraints  need  to  be  imposed  to  obtain  a  unique 
solut ion. 

As  early  as  1954  Bracewell  and  Roberts (ref. 1)  proposed  an  iterative 
deconvolution  technique  which  progressively  sharpens  the  estimated  distribution 
and  ultimately  converges  to  a  unique  solution  -  termed  the  'principal  solution'. 
This  solution  has  the  additional  property  that  it  contains  no  wave-number 
components  which  give  a  zero  output  when  convolved  with  the  array  response. 
Although  this  'principal  solution'  has  greater  resolution  than  the  conventional 
beamformer  estimates,  its  resolution  is  still  limited  and  it  also  suffers  from 
ringing,  i.e.  an  increase  in  the  sidelobe  levels.  Its  effect  in  terms  of  the 
spatial  autocorrelation  function  can  readily  be  seen.  For  a  continuous  line 
array  the  method  replaces  the  Bartlett  (or  triangular)  window,  multiplying  the 
spatial  autocorrelation  distribution  by  the  square  window.  The  effect  of  this 
is  well  known;  the  resolution  is  doubled  at  the  expense  of  increased  sidelobe 
levels. 

However,  the  technique  has  received  considerable  attention.  Axelrod  et  al 
(ref.l)  have  shown  how  the  method  can  alternatively  be  formulated  in  terms  of  a 
least  squares  expansion  of  the  wave -number  power  spectrum  using  a  set  of 
Fourier  coefficients.  Anderson  and  Tittle(ref . 1) ,  by  defining  the  principal 
solution  mathematically  have  shown  that  the  principal  solution  only  retains  a 
finite  numbeT  of  terms  in  the  expansion.  Furthermore  all  these  terms  corres¬ 
pond  to  components  which  oscillate  in  wave-number  (i.e.  essentially  angle) 
below  a  certain  frequency. 

Further  work  by  Wilson(ref . 1)  and  in  particular  McDonough(ref . 1)  has 
developed  the  method  and  applied  it  to  some  numerical  examples.  McDonough  has 
also  generalised  the  iterative  method  of  Bracewell  and  Roberts  to  an  array  of 
arbitrary  geometry. 

The  techniques  discussed  all  suffer  from  one  common  limitation;  i.e.  the 
deconvolution  is  effected  using  the  beam  powers.  Thus  the  relative  phase 
information  of  either  beams  or  receiver  outputs  is  lost.  A  method  overcoming 
this  limitation  has  been  proposed  by  Nuttall(ref .1)  whereby  the  incident  field 
(i.e.  the  noise  field),  in  a  similar  manner  to  Axelrod  et  al,  is  expanded  in  a 
finite  set  of  Fourier  coefficients.  However,  instead  of  forming  a  least 
squares  approximation  using  beam  powers,  Nuttall's  method  uses  the  crosspower 
spectral  matrix  of  the  receiver  outputs.  Yen (ref. 2)  has  also  adopted  this 
approach,  applied  it  to  a  linear  array  and  obtained  an  expression  for  the 
coefficients  of  the  Fourier  expansion  in  terms  of  the  relative  phase  delays  of 
the  receivers. 

The  expansion  of  the  noise  field  in  a  Fourier  series  in  wave-number 

(i.e.  y  cos  also  deserves  some  comments.  In  this  approach  the  lowest  order 

term  represents  the  isotropic  noise  component  and  to  represent  anisotropic  noise 
fields  higher  order  Fourier  terms  are  required.  In  the  limit  a  single  plane 
wave,  represented  by  a  delta  function  distribution  of  the  noise  field,  requires 
an  infinite  number  of  coefficients  for  exact  representation.  The  basis  of  this 
paper  is  to  adopt  what  can  be  conceptually  thought  of  as  being  almost  the  exact 
opposite  to  the  above  approach,  that  is,  the  incident  field  at  the  frequency  of 


interest  is  assumed  to  be  composed  of  N  independent  plane  waves  (i.c.  N  delta 
functions)  from  predetermined  directions  and  the  field  is  assumed  homogeneous. 
This  reduces  the  estimation  problem  to  one  of  estimating  only  the  complex  beam 
amplitudes  or  powers  using  an  array  of  K  discrete  receivers. 

This  assumption  has  been  successfully  exploited  by  d'Assumpcao(ref .2)  to 
derive  some  new  quadratic  estimators  based  on  maximum  likelihood  and  minimum 
variance  criteria. 

Under  these  assumptions  the  phase  delays  corresponding  to  the  selected 
directions  may  be  incorporated  in  the  deconvolution  process  which  is  represented 
in  matrix  formulation  since  both  the  number  of  receivers  and  arrivals  is  finite. 
The  method,  although  formulated  for  a  two-dimensional  noise  field  may  readily  be 
extended  to  three  dimensions.  In  Section  2  an  iterative  deconvolution 
technique,  which  is  a  matrix  extension  of  the  method  used  by  Bracewell  and 
Roberts,  is  formulated  under  the  above  assumptions.  Two  similar  techniques  are 
proposed  which  use  the  array  response  at  the  frequency  of  interest  to  deconvolve 
either  the  complex  beam  outputs  or  the  beam  powers.  The  assumption  that  the 
directions  of  the  N  plane  waves  are  known  enables  the  total  leakage  in  any  beam 
to  be  estimated  and  subtracted  out.  The  iterative  method  is  used  to 
successively  refine  and  subtract  out  the  leakage  in  all  beams  from  the  other 
beams.  It  is  then  shown  in  Section  3  that  the  limits  of  these  iterative  methods 
converge  to  a  generalization  of  some  particular  quadratic  estimators  derived  in 
reference  2 

Recently,  Yen(ref.3)  using  a  similar  assumption  for  the  incident  noise  field 
has  used  the  Prony  method  to  estimate  the  direction  and  powers  of  the  incident 
N  plane  waves  for  a  line  array.  The  expressions  derived  in  this  report 
together  with  the  examples  of  Section  5  can  be  used  to  show  an  equivalence  of 
the  power  deconvolution  method  to  the  linearized  part  of  the  Prony  method  used 
by  Yen. 

In  Section  4  the  suitability  of  these  estimators  in  the  presence  of  noise 
with  an  unknown  covariance  matrix  is  discussed  and  some  general  results  are 
proved. 

Finally  some  examples  of  the  application  of  the  method  to  a  line  array  of  K 
equi spaced  receivers  are  given. 


2.  FORMULATION  OF  THE  ITERATIVE  METHOD 


he  assumption  that  the  incident  field  consists  of  N  independent  plane  waves 
enables  x^ ,  j  =  1,  2,  . . . ,  K  (the  complex  spectral  amplitude  of  the  output  of 

the  j-th  receiver  at  a  frequency  f)  to  be  expressed  as  a  linear  combination  of 
y^  (the  complex  spectral  amplitude  of  the  k-th  plane  wave  signal  at  some 

arbitrary  reference  point) .  The  physical  geometry  of  the  array  determines  this 
particular  combination.  It  follows  that: 


N 


k=l 


where  =  exp(i^^)  and  0^  are  the  phase  delays  corresponding  to  a  signal 

from  the  k-th  direction  at  the  j-th  receiver.  For  the  example  of  a  line 

sin  0. 


array  of  equispaced  receivers  0^  *  2ff  fj 


d  where  d  is  the  separation 
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of  adjacent  receivers.  Denoting* 


X  -  (Xj i  X^,  . . . |  xKj  , 


y  =  C/j,  y2»  yN)  , 


V  -  {_  Vj 


jk  ' 


a  K  x  N  complex  matrix,  equation  (1)  can  be  rewritten  as: 

x  =  Vy  .  ( 

2.1  Deconvolution  of  complex  beam  outputs 

The  output  of  a  conventional  beamformer,  denoted  as  y^  is  then  given 


v(0)  _  _A 

y  NK 


In  order  to  deconvolve  the  yj  },s  the  method  of  progressive  substitutions 

as  discussed  by  Bracewell  and  Roberts  is  used.  However,  since  the  x^'s  are 

assumed  related  to  the  y^'s  by  equation  (2)  it  is  now  possible  to  deconvolve 

not  just  with  respect  to  the  array  window  but  also  with  respect  to  the  phase 
delays  <t> .  . . 

13  (01 

If  the  initial  estimate,  yv  J  ,  was  the  incident  noise  field  distribution 
then  the  receiver  outputs  would  be 


x  -  Vyw  . 

Hence  the  output  of  the  conventional  beamformer  would  be  Wy^  where 
V^V 

W  =  -jjjj—  .  A  measure  of  the  error  in  the  original  estimate  y ^  is 

then  defined  by: 


«<°>  =  y<°>  -  ny«°) 


*  T  denotes  the  transpose  and  H  the  complex  transpose  of  either  a  vector  or 
a  matrix. 

*  The  choice  of  NK  rather  than  K  as  the  normalizing  factor  ensures  convergence 
of  the  iterative  series  (see  Appendix  II)  at  the  expense  of  biasing  the 
power  spectral  estimates.  The  reduction  of  this  bias  is  discussed  in 
Section  3. 
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A  new  approximation  to  the  incident  distribution  is  obtained  by  correcting 
for  this  error,  viz: 

yt«  =  y««  . 

In  general,  let  y^  represent  sane  n**1  order  approximation  to  the 
incident  distribution.  If  this  represented  the  true  distribution  then  the 

output  of  a  conventional  beamformer  would  be  lVy^n^.  The  error  term  between 

this  output  and  the  observed  output  of  the  beamformer,  y^,  is  denoted  as 

e ^  and  is  defined  by: 


e  =  yW)  -  Wy™  . 


(3) 


The  method  of  progressive  substitutions  then  implies  that  a  new 
approximation,  y^n+^,  to  the  true  incident  distribution  can  be  chosen  as: 


y(n+l)  ,  y(n)  +  e Cn)  # 

This  reduces  to: 

y<"*»  •  y<°)  .  (I  -  W)y'n> 


(1) 


by  virtue  of  equation  (3) . 

An  alternative  formulation  of  the  iteration  which  gives  the  same 
results  but  has  a  direct  physical  significance  will  now  be  given. 

Consider  the  ij**1  element  of  W  for  i  =£  j ;  it  represents  the  (biased) 

response  of  the  i  beam  to  a  plane  wave  from  the  j  direction.  Thus  a 
given  row  of  NW  is  the  amplitude  polar  response  of  the  array  evaluated  at 
the  wave-numbers  (or  angles)  corresponding  to  j  =  0,  1,  .  •  •  N-l.  Now 

suppose  y^  is  some  approximation  to  the  incident  distribution.  The 

components  which  distort,  through  leakage,  an  estimate  of  i^  beam  will  be 
given  by: 


W.  .y 
irj 


(n) 


This,  for  all  beams,  reduces  to 


(W  -  I)y^ 


since,  W. .  =  1/N.  In  order  to  attempt  to  cancel  the  effect  of  these 

11  fO) 

sid Globes  the  'leaked'  components  can  be  subtracted  from  yv  1  ;  the 

conventional  estimate  of  the  spectrum.  This  leads  to  a  better  approxima¬ 
tion  of  the  wave-number  spectrum  denoted  as  y*-n+^  which  is  given 
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by  equation  (3).  The  convergence  rate  of  the  iteration  will  be  increased 
by  ensuring  that  the  off-diagonal  elements  of  W  are  small.  This  is 
alternatively  interpreted  as  requiring  the  array  polar  diagram  to  have 
nairow  beamwidth  and  low  sidelobes. 

2. 2  Deconvolution  of  beam  powers 

In  this  section  the  iterative  technique  is  used  to  deconvolve  the  beam 
powers.  The  conventional  estimate  of  the  N  beam  powers  s ^  is  defined  by: 


<  (VHxx**V)  .  .> 

_ _  li 

Ink!2 


(5) 


where  <> denotes  the  ensemble  average  (in  practice  replaced  by  a  time 
average) . 

Define  the  vector  m  by: 


m  =  <  x  ®  x*> 

where  ®  denotes  the  direct  products  and 

A  =  V  O  ^  (6) 

where  ©  is  the  Khatri-Rao  product  (i.e.  (A  ©  B)  _  =  a^b^  ^  where 

i  =  (k  -  1)K  +  l )  and  the  dimension  of  A  is  KJ  x  N.  Then  equation  (5) 

for  s.  can  be  written  as: 
i 

/n(O)  =  AHm 

Ink!2 

Also  defining: 

S  =  <yyH> 


and  substituting  in  equation  (5)  for  x  it  follows  that: 

(v^Vsv^V) . . 

A  _  _ _  11 

Si  |  NKj 2 

=  (WSW)..  . 


Furthermore  the  assumption  that  the  signals  from  the  differing  directions 
are  uncorrelated  implies  that  S  is  diagonal  and  so  the  above  equation 
reduces  to: 


s 


(0) 


W  □  WTs 


where* 


and 


s. 


1 


=  S. . 
11 


(W  □  WT).. 


CV^yiyVy 

[nkP 


(7) 


Thus  an  initial  approximation  to  s  can  be  s ^  which  would  result  in  an 
error  term  of: 


e(°>  .  <8> 

which  can  be  used  to  correct  the  original  estimate,  i.e. 

s(1)  =  s(0)  ♦  .<«  . 

This  process  is  now  repeated  and  the  general  iterative  equation  becomes: 
s(n+1)  =  s(0)  +  (i  .  W  □  WT)s^n^  . 

T 

As  in  the  previous  section  the  iteration  matrix,  i.e.  W  E3W  ,  has  a 

T 

physical  interpretation.  From  equation  (7)  it  follows  that  N2W  G)W  is 
simply  the  polar  diagram  of  the  beam  powers  and: 


)  (W  GJ  WT)  .  .sfn) 

ij  J 

i.L 

is  the  leakage  into  the  i  beam  of  the  powers  from  the  N-l  other 
directions.  Thus  the  technique  may  be  considered  as  reducing  leakage  of 
powers  (either  because  of  a  broad  beamwidth  or  high  side  lobe  levels)  from 
one  beam  to  another  by  using  the  'a  priori'  knowledge  of  the  array's  polar 
response. 


*  The  Hadamard  product,  □  ,  of  two  matrices  A  and  B,  is  defined  by: 


3.  SOLUTION  AND  CONVERGENCE  OF  THE  ITERATIVE  METHOD 


In  this  section  and  in  the  remainder  of  this  paper,  let  U  denote  either  of 

the  N  x  N  matrices,  i.e.  W  or  II  Q  and  let  z^  represent  either  y^  or  s^. 
From  either  equation  (4)  or  (7)  it  follows  that: 

z(n)  =  (I  +  (I  -  U)  +  ...  (I  -  U)n)z(0^  (9) 

and  (see  Appendix  I)  this  can  be  shown  to  reduce  to: 

z(n)  =  (1  -  (I  -  U)n+1)U"z(0)  , 


where  U  is  any  generalized  inverse (ref .4)  of  the  matrix  U.  Unfortunately 
convergence  of  this  series  does  not  hold  in  general  since  the  eigenvalues  of 
(I  -  U)  are  not  all  less  than  1. 

3.1  Convergence 

A  simple  modification  of  the  recurrence  relation  allows  convergence  of 
the  iteration.  Replace  equation  (9)  by: 

,'■1  .  .  CXI  -  lOz'"-1’  . 

The  effect  of  the  A  is  to  modify  the  recurrence  formula  to: 

z(n)  =  Az^  ♦  e<n"1> 


where  e ^  is  the  error  term  defined  by  either  equation  (3)  or 
equation  (8).  As  shown  in  Appendix  B  the  restriction  0  <X  <  1  will 
guarantee  convergence  of  the  series.  This  requirement  also  has  a 
heuristic  physical  interpretation  since  the  condition  that  A  <  1  can  be 
considered  as  automatically  allowing  for  the  fact  that  for  N  >  K  there 
must  always  be  leakage  from  one  beam  to  another.  This  is  a  reflection 
of  the  fact  that  it  is  impossible,  for  N  >  K,  to  steer  more  than  K-l  nulls. 
As  a  special  case,  for  N  =  K  and  V  non-singular,  convergence  of  the  series 
for  A  =  1  is  possible  since  it  is  now  possible  to  steer  K-l  nulls  and  so 
prevent  any  leakage  into  a  selected  direction. 

3.2  Limit  (beam  amplitudes) 

The  matrix  identity 


(V(Al  -  U)^I+(Xl  -  U)  +  (Al  -  U)2+  ..  +  (Xl  -  U)n^  =  (V(XI  -U)"+i)  (10) 

always  holds.  Since  A  <  1  and  W  =  (V^V)/(KN)  it  follows  that  (l-A)l+w 
is  always  non- singular.  It  then  follows  that: 

=  j  I  -  (XI  -  10"*1]  \  (1  -  X)l  ♦  I*  I'M01  ■ 
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However,  since  the  eigenvalues  of  Xl  -  W  are  always  strictly  less  than  1 
(see  Appendix  B)  it  follows  that: 


lira  (Xl  -  W) 
n+6° 


,  •  (n) 

lim  yK  ’ 


(1  -  X)i  ♦  w  “ 1y(-0) 


v”v  1  -1  VHx 


which,  in  general  will  give  different  estimates  for  differing  values  of  X. 
However,  from  reference  3,  the  Moore-Penrose  pseudoinverse  of  V,  denoted 


as  V+  is  defined  as: 


V+  =  lim  (5 1  +  . 

5*0 


It  then  follows  that: 


lim  lim  y 
X*  1  n*°° 


(n)  _ 


=  V  x  . 


Thus,  as  a  generalization  of  conventional  beamforming  V**  is  replaced  by 
V+  the  Moore-Penrose  pseudoinverse. 

It  is  worth  noting  that  any  choice  of  scaling  for  y^  which  ensures 
convergence  gives  rise  to  the  same  limiting  solution  V  x.  However  as 
discussed  in  Section  2  the  above  solution  is  biased.  A  sensible  con¬ 
straint  to  further  impose  is  that  all  beams  have  the  same  maximum  response 
in  their  look  directions.  Since  the  output  of  the  inverse  processor  in 

the  n**1  direction  due  to  a  plane  wave  incident  from  that  direction  is: 

<V+Vn  * 


where  vn  is  the  n**1  column  of  V,  it  follows  that  for  unity  response  in 
the  look  direction  the  appropriate  form  is: 


(V  x). 


C v+x)i 

yi  =  ~ — 
(V  V).. 


3.3  Limit  (beam  powers) 

Since  the  matrix  identity  (equation  (10))  always  holds  it  follows  that: 


s(n)  =  -  (Xi  -  w  □  wH)n+1j  j^(i  -  X)i  +  w  □  iv9  J  1s(0) 

provided  the  inverse  of  {(1  -  X)I  +  W  □  W**]  exists.  The  further  matrix 
identity(ref . 5)  : 


<VHV)EKVTV*)  =  (V  o  V*)H(V  G  V*)  (11) 

guarantees  that  an  inverse  exists  for  X  <  1.  Consequently  since  the 
eigen-values  of  (Xi  -  W)  are  always  strictly  less  than  1  it  follows  that: 

lim  (Xi  -  W  □  WT)n+1  =  0 

n+°° 

and  so 


lim  s^  =  (1  -  X)I  -  W  □  W 

n>°°  L 

Using  equation  (11)  this  reduces  to 


5  s<°> 


lim  s(n) 

TT»°° 


(1  -  X)I 


A 

IKN}1 


-1 


Am 

\m* 


where  A  =  V  G  V*.  As  in  the  previous  case  the  limiting  form  as  X-*-l“  is 
the  pseudoinverse  and  so  an  'optimum'  processor  is: 


s 


A  m 


and  substituting  for  A  and  m  this  becomes: 

s  =  (V  O  V*)+  <x  8  x*> 

As  before  the  limiting  solution  is  independent  of  the  scaling  of  s ^  and 
a  solution  which  allows  unity  response  in  the  look  directions  is  given  by: 


(A+m). 

(A+A) . . 
13 
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4.  DISCUSSION 

4.1  Complex  beam  outputs 

The  estimator  y  =  V+x,  termed  the  inverse  estimator,  has  a  number  of 
interesting  properties  which  relate  it  to  some  quadratic  processors  which 
have  recently  been  proposed. 

(a)  If  V  is  square  and  non-singular  (i.e.  K  =  N)  then  V+  =  V1.  For  a 

chosen  direction  the  processor  V  steers  nulls  in  the  remaining 
(K  -  1)  directions  and  thus  forms  an  unbiased  estimator  of  the  wave- 
number  spectrum.  This  estimator  has  been  derived  in  reference  2 
using  a  maximum  likelihood  technique. 

(b)  If  for  an  arbitrary  N,  the  K  rows  of  V  are  linearly  independent 

then  V+  =  V^(W^)  *.  Apart  from  the  scaling  factor,  (V+V)  ,  this 

is  the  minimum  bias  estimator  derived  in  reference  2.  nn 

(c)  In  general  a  class  of  solutions  of  the  equation: 

x  =  Vy  (12) 

where  W  y  =  y,  is  given  by: 

y  =  V'x 

where  V  is  any  generalized  inverse.  The  solution  derived  here, 

V+,  has  the  additional  property  that  of  all  the  y's  which  satisfy 
equation  (12),  it  is  the  one  with  the  minimum  norm(ref.4).  That 

is  llyll2  =  S  |  y^| 2  is  also  minimized.  Alternatively  if  the 
i 

receiver  outputs,  i.e.  the  x^,  are  uncorrelated  then  the  pseudo¬ 
inverse  is  the  one  which  also  minimizes  the  norm  of  the  weight 
vectors  and  thus  limits  the  superdirectivity  of  the  array. 

(d)  In  general  the  effect  of  noise  (and  also  signals  arriving  from 
directions  not  accounted  for  by  the  v.,  's)  will  be  such  x  is  only 

approximately  equal  to  Vy.  The  choice  of  y  =  V  x  as  an  estimator 

is  still  an  appropriate  one  since  it  is  shown  in  reference  4  that 
it  is  the  vector  of  minimum  norm  which  minimizes  llx  -  Vyll2  .  That 
is,  the  estimator  y  =  V  x  is  the  best  (in  a  least  squares  sense) 
plane  wave  solution  to  the  problem  subject  to  the  constraints 
imposed  by  the  assumed  source  directions. 

U 

(e)  The  total  power  output  from  the  receivers  is  x  x  .  If  y  is  defined 

as  y  *  Ux  then,  when  U  is  any  unitary  matrix  (lAj  =  I),  the  total 

power  in  the  y's,  yHy,  equals  that  in  the  x's.  This  is  analogous 
to  Parseval's  theorem  in  Fourier  analysis  and  is  an  expression  of 
the  conservation  of  energy.  If  the  V's  are  not  unitary  then  this 
equality  does  not  hold  but  a  sensible  criterion  would  be  to  require 
that  it  holds  as  an  approximation.  This  then  implies  that: 

lixHx  -  yHy)l2 


should  be  minimized. 


(13) 


11 
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and 

(V+V)H  =  V+V  . 

Differentiating  the  above  with  respect  to  z  and  equating  to  zero 
then  implies  that  equation  (14)  is  minimized  when: 

z  =  V+Vz  . 

Hence  the  expression  for  y  reduces  to: 

y  =  V  x  . 

Thus  the  Moore-Penrose  pseudoinverse  is  that  solution  minimizing: 

it  H  aHah2 
llx  x  -  y  ylr 


i.e.,  it  most  closely  conserves  the  power. 

An  important  consequence  of  this  result  is  that  if  there  is  a 
strong  source  in  the  x’s  with  a  wave-number  not  accounted  for  by 
the  V's  then  there  will  be  considerable  leakage  of  the  power  of 
this  source  into  the  estimated  y^.  The  (dis) advantage  of  this 

will  depend  on  whether  this  strong  source  is  either  a  desired 
signal  or  an  interference. 

(f)  Idempotency  of  solution 

Suppose  the  inverse  estimates,  i.e.  y  =  V+x,  are  the  true 
distribution  of  sources.  The  receiver  outputs,  x,  become: 

x  *  VV+x 


* 
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and  the  output  of  the  inverse  estimator  is: 

V*x  =  V+W+x 
=  V+x  . 

This  property  of  the  inverse  estimator  is  termed  idempotency 
and  it  is  precisely  the  lack  of  this  property  for  the  conventional 
processor  that  allows  the  iteration  to  be  effected. 

4.2  Beam  powers 

(a)  If  V  is  square  and  non-singular  (i.e.  K  =  N)  and  since: 

w  □  WT  -  (V  G  V*)H(V  o  V*) 

Ink!2 

T 

then  W  □  W  is  positive  definite  and  hence  non-singular.  Thus 
the  series  is  convergent  for  X  =  1  and  the  limit  becomes: 

(W  □wT)"1s(0) 

which  on  substitution  of  equations  (6)  and  (11)  becomes: 

jNKj2  (AHA)-1s(0) 

which  has  been  derived  in  reference  3  as  a  least  squares  estimator 
(see  also  below)  under  some  more  general  conditions. 

(b)  From  Section  3  the  approximation: 

m  *  VSV9 

where  S  is  diagonal  can  be  rewritten  as 

m  ^  As 

where  s  is  the  vector  of  the  diagonal  elements  of  S.  As  in  the 
previous  section  a  particularly  suitable  approximation  is: 

s  =  A  m 

since  in  addition  to  minimizing: 

II  m  -  Asll2  (or  dm  -  VSVH|lJ ) 


13  - 
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it  also  minimizes: 


and  hence  limits  the  superdirectivity  of  the  estimator. 


5.  EXAMPLES 

Consider  a  linear  array  of  K  equispaced  receivers  separated  by  a  distance  d. 
If  the  incident  distribution  is  assumed  to  be  two-dimensional  and  composed  of 
N  plane  waves  with  wave  numbers  k.  (=  (2it  sin  0  ,)/X)  then  the  matrix  of  phase 
delays,  V,  has  the  form:  ■* 


V..  =  zt 

il  J 


where 


2flidk. 


z.  -  e 


This  is  a  Vandermonde  matrix(ref . 6)  and  us 2  may  be  made  of  its  special 
properties.  Furthermore,  if  the  k^  are  chosen  to  correspond  to  N  arrival 

directions,  equispaced  in  sin  0,  and  lying  between  ±w/2  then  V  is  simplified  to: 


V..  =  z 

ij 


where 


z  =  e 


4*  id  An 


For  such  a  V,  W,  defined  as  (V^V)/(KN)  and  the  matrix  to  be  used  in  the 
amplitude  iteration,  is  given  by: 

1  .  1) 


w  -  1) 


where  u  =  p  -  p.  The  matrix  W  □  W1  can  be  readily  reduced  to: 


sin2  ^  (p  -  P) 


I  sin  ^  IP  -  J 
sin*  0>  ■  » 
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which,  apart  from  the  factor,  is  the  polar  response  of  a  line  array  of  K 

equispaced  receivers.  Since  both  W  and  W  □  W  are  Toeplitz  it  then  follows 
that  the  term: 


U  z 

pv  v 


in  the  iteration  becomes: 


Wp 


This  is  simply  the  result  of  the  first  N  terms  of  the  cyclic  convolution  of  the 
vector: 


(u0>  uj» 


»  uM-l*  U-N+1*  u-N+2’ 


• t  $  0) 


with  the  vector: 


( z(n)  z(n) 

Vzo  »  Z1 


'  ft  ft 

•»  ZN-1'  u'  u 


As  a  result,  if  N  is  chosen  to  be  a  power  of  two,  then  the  convolution  may 
be  efficiently  evaluated  by  use  of  the  fast  Fourier  transform.  Unfortunately 

the  zjn+1)  for  j  =  N,  ...,  2N-1  produced  by  this  cyclic  convolution  are  not 

zero  and  so  the  iteration  cannot  be  effected  completely  in  the  transform  space. 
This  can  alternatively  be  realized  by  observing  that  U*  is  not  in  general, 
Toeplitz.  (Aside:  the  method  proposed  by  Bracewell  and  Roberts  assumed  that 
the  iteration  could  be  effected  completely  in  the  transform  domain.  This 
amounts  to  assuming  the  array  distribution  to  be  cyclic  and  so  is  why  they  only 
succeeded  in  deconvolving  the  array  shading) . 

V  is  generally  of  full  rank  and  so  V  x,  the  limit  of  the  amplitude 
iteration,  may  be  written  as: 


VH(WH)'1x  . 


In  particular: 


cw\„  -  :  R 


and  is  a  complex  Toeplitz  matrix. 
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Efficient  algorithms  exist  for  the  inversion  of  this  matrix  enabling  V+  to 

be  rapidly  calculated.  For  the  particular  case  of  K  =  N  then  V+  =  V~  and 
once  again,  since  V  is  a  Vandermonde  matrix,  it  may  be  efficiently  inverted. 

H 

Similarly  A  A  can  be  shown  to  be  a  Toeplitz  matrix  and  is  non-negative  since 
all  its  elements  are  positive. 

Some  examples  of  these  for  various  values  of  N,  K  and  d/X  will  be  given. 

5.1  K  receivers  and  1  plane  wave 

In  this  case  V  is  the  K  x  1  column  vector  and  V+  trivially  reduces  to 

v“ 

-jr  which  is  the  conventional  processor.  The  iterative  equation  for  the 
amplitudes  reduces  to 


y(n)  =  y(0)  +  (X  -  i)y 


(n-1) 


which  for  X  =  1  becomes  y 


(")  =  y(0) 


as  would  be  expected. 


5.2  K  receivers  and  N  sources  at  half  wavelength  spacing 

When  d  »  X/2  it  follows  that  z  =1  with  the  result  that 


VV*1  = 


NI,. 


where  IK  is  the  K  x  K  identity  matrix.  As  a  consequence  the  first  step 
in  the  iteration  equation: 


y'°>  .  (xi 


vV\  A 

NK  /  NIC 


reduces  to: 


For  X  =  ^  the  first,  and  consequently  all  successive  iterations,  reduce  to 

the  conventional  beamformer.  It  also  follows  that,  for  X  =  1,  either  as  a 

4  yH 

consequence  of  the  iteration  or  the  fact  that  V  =  —  ,  the  limit  y  reduces 

VHx  N+  K 

to  the  biased  estimator,  -y  .  Furthermore  since  (V  V)^  =  -yj-  the 

A 

standard  solution  y  ,  corresponding  to  unity  response  in  the  look 
direction,  is  obtained. 

5.3  K  receivers  and  K  sources  at  quarter  wavelength  spacing 

n  if 

The  number  of  independent  beams  lying  between  -  and  —  for  a 

conventional  processor  at  a  quarter  wavelength  is  y  .  The  first  of  these 

>1 

is  plotted  in  heavy  lines  in  figure  1.  The  extra  y  beams  incorporated 
in  V  are  redundant  beams  spaced  half  way  between  the  adjacent  independent 
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beams;  the  first  redundant  beam  is  plotted  in  dashed  lines  in  figure  1. 
Now  w  (=W  )  ,  essentially  the  array  polar  diagram  is  given  by: 

r  r*  +r 


for 

P 

odd 

P  N1  (zP  -  1) 

=  0 

for 

P 

even  =£  0 

J_ 

M 

for 

P 

=  0. 

Consider  the  leakage  of  y^  into  the  other  beams; 


from  figure  1  zero  power 


will  leak  into  the  even-number  beams  since  they  correspond  to  nulls  in  the 
polar  diagram  (this  is  reflected  in  the  fact  that  w^  =  0  for  p  even).  The 

power  leakage  into  the  odd  beams  is  given  by  w  y  where  w  is  given  by  above 

and  y^  is  the  output  of  the  zero1-  beam.  Thus  the  form  of  w^  enables  the 


leakage  from  y^  (and  in  general  all  the  beams)  to  be  removed  from  the  other 
beams  as  discussed  in  Section  2. 

Some  other  points  follow  from  this  interpretation  of  the  deconvolution 
technique.  In  figure  1  the  side  lobes  of  any  beam  are  shown  extending  into 
what  is  termed  the  'non-physical  region'.  This  region  corresponds  to  plane 
wave  disturbances  either  generated  in  the  array  of  receivers  ot  propagating 
across  the  array  with  a  velocity  less  the  cne  assumed  in  calculating  the 
phase  delays.  It  is  an  inherent  assumption  of  these  examples  that  the 
amplitude  of  any  of  these  effects  is  zero  or  at  least  very  small.  Some  of 
these  wave-number  beams  could  be  incorporated  in  V  but  unfortunately  any 

N 

attempt  to  account  for  the  full  —  (at  a  quarter  wavelength)  independent  ones 

would,  since  the  matrix  W  becomes  circular,  reduce  V+  to  V^;  the  conven¬ 
tional  processor.  Once  again  this  can  be  intuitively  seen  from  a 
consideration  of  the  zeros  of  the  polar  diagram. 


6 .  CONCLUSIONS 

Estimates  of  the  angular  distribution  of  the  power  incident  on  an  array  made 
by  using  the  outputs  of  a  conventional  frequency-domain  beamformer  are  distorted 
by  leakage  from  one  beam  to  another.  Two  techniques  have  been  proposed  in  this 
paper  which  use  a  knowledge  of  the  polar  response  of  the  array  to  minimise  this 
leakage.  The  techniques  are  conceptually  similar,  the  difference  being  that 
one  uses  the  complex  beam  outputs  whereas  the  other  uses  the  beam  powers. 

Both  techniques  can  be  effected  as  a  series  of  iterative  deconvolutions  using 
either  the  narrowband  array  amplitude  or  polar  response  and  the  narrowband  beam 
outputs  or  powers  respectively  of  a  conventional  beamformer.  Alternatively 
either  the  receiver  outputs  or  the  receiver  crosspower  spectral  matrix  may  be 
used  directly  to  evaluate  the  limits  of  the  two  iterative  methods.  In  practice 
the  choice  would  be  determined  by  implementation  requirements. 

A  particular  attraction  is  that  the  techniques  are  linear  and  so  the  concept 
of  a  polar  diagram  is  useful.  The  convergence  of  the  iteration  as  a  function 
of  the  parameter  X  and  the  stability  of  the  matrices  to  be  inverted  are  areas 
warranting  further  investigation. 
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APPENDIX  I 


Given: 


A  VECTOR  IDENTITY 


z 


(0) 


Uz 


Cl.  ID 


A  (°) 

and  z 


is  defined  by  the  recursive  relation: 


zW  -  z<°>  .  (I  -  Ujz'"-1' 

it  is  required  to  show  that: 

z(n)  =  (I  -  (I  -  U)n+1)lTz^ 


(1.2) 


(1.3) 


where  U  is  any  generalized  inverse  of  U. 

Before  proceeding  with  an  inductive  proof,  it  follows  directly  from 
equation  (1.1)  and  the  property  of  UU  U  =  U  for  any  generalized  inverse  that: 


UlTzW  =  *«» 


(1.4) 


where 


z(°)  .  y(0)  or  s(0) 


Assuming  that  (1.3)  holds  for  z^n  ^  and  then  substituting  for  z^n , 
equation  (1.2)  reduces  to: 


z(n)  =  z(0)  +  (i  .  u)(I  -  (I  -  U)n)U'z(0) 

=  z<°>  -  UU'z(0)  +  (I  -  (I  -  U)n+1)lTzW 
=  (I  -  (I  -  U)n+1)U'z(0) 


where  the  last  step  follows  directly  from  (1.4).  To  complete  the  induction 
it  is  necessary  to  prove  equation  (1.3)  when  n  =  0.  For  n  =  0  equation  (1.3) 
reduces  to: 


z‘°>  =  UU-z<°> 

which  follows  directly  from  (1.4).  .  . 

The  proof  is  completed  by  showing  that  z^  defined  by  equation  (1.3)  is 
unique  for  any  choice  of  U”. 

Now  any  generalized  inverse  of  U,  Ui  ,  can  be  written  as 
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(o) 

(o) 


,(o) 


Hence  equation  (1.3)  is  unique  if 


(I  -  (I-U)  j  (I  -  U'U)  =  0  . 


Bur 


(I  -  (I-U)n)  (I  -  U'U)  .  (I  +  (I-U)+  ...  +  (I  -  U)n-1)U  (I  -  U'U) 
Thus  defined  by  equation  (1.3)  is  unique. 
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APPENDIX  II 

CONVERGENCE  OF  A  SERIES 

In  order  to  show  that  the  geometric  series  I  +  (Xl  -IV)  +  ...  +  (Xl  -  W)n 

converges  as  n  +  00  it  is  necessary  to  show  that  lim  (Xl  -  W)n  as  n  *  00  is  zero. 
A  necessary  and  sufficient  condition  for  this  is  that  all  the  eigenvalues  of 
Xl  -  W  are  less  than  1. 

v^V 

Let  P(<  K)  be  the  rank  of  W  a  K  x  K  matrix  defined  by  W  =  where  V  is 

any  K  x  K  matrix.  It  always  holds(ref.5)  that: 

W  =  QH  A  q  (II.  1) 


where 

A  = 

diag  |  X^,  \22,  . . . ,  Xp,  0,  0,  . . . ,  0 

and 

Q  = 

is  a  unitary  N  x  N  matrix. 

Convergence  of  the  geometric  series  I  +  (Xl  -  W)  +  ...  +  (Xl  -  W)n  is 
guaranteed  provided: 

IX  -  X.l  <  1  i  =  1,  2,  ....  P 


and 


I  X|  <  1 


These  two  conditions  can  be  seen  to  be  satisfied  provided  0  <  X  <  1  and 
|XJ<  1.  The  first  condition  can  easily  be  satisfied  by  an  appropriate  choice 

of  X. 

If  the  eigenvalues  of  the  matrix  V*V  are  it  follows  that  the  eigenvalues 

of  L  '/V  are  less  than  or  equal  to  unity.  Since  yHy  is  positive  definite, 
'max 


i.e.  M.  =  X.  it  follows  that: 
1  1 


max 


r 

E 


i=l 


However,  it  holds  that: 


Mi  =  Tr(\,HV) 


i*l 
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